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ABSTRACT 

A  computer  program  was  developed  for  solving  equality 
and  inequality  constrained  optimization  problems  by  the 
Augmented  Lagrange  Multiplier  method.   The  program  was 
developed  specifically  for  use  in  engineering  design. 

The  historical  evolution  and  theoretical  development  of 
the  multiplier  method  is  presented.   Several  examples  are 
used  to  demonstrate  the  effects  of  penalty  parameters  and 
multipliers  on  the  convergence  and  accuracy  of  the  method. 
Computational  experience  with  variations  to  the  method  is 
documented. 

A  brief  literature  search  of  the  multiplier  method's 
application  to  engineering  design  is  summarized.   The  method 
is  demonstrated  with  several  mathematical  and  engineering 
examples.   A  comparison  to  classical  penalty  methods  and  the 
method  of  feasible  directions  was  performed  in  each  case. 
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I.   INTRODUCTION 

Several  methods  currently  exist  which  solve  the  function 
minimization  problem.   While  each  has  its  advantage  in  special- 
ized cases,  no  one  method,  including  the  multiplier  method,  is 
universally  more  efficient  and  accurate  than  the  others.   The 
Augmented  Lagrange  Multiplier  method  has  been  shown  to  be  an 
efficient  and  accurate  method  of  engineering  design  optimiza- 
tion.  The  multiplier  method  is  particularly  advantageous  in 
nonlinear  equality  constrained  problems.   It  has  become 
increasingly  popular  in  engineering  design  optimization 
because  of  its  good  convergence  rate  and  its  attractive 
theoretical  properties. 

The  general  nonlinear  problem  is  defined: 

Minimize  f(X)  (1.1a) 

subject  to  g.  (X)  <  0  (i=l,...,J£)  (1.1b) 

and  h.(X)  =  0  ( j=l , . . . ,m<n)  (1.1c) 

and  X£  <  X.  <  X^  (i=l,...,n)  (l.ld) 

where  X  is  a  vector  of  n  design  variables,  f (X)  is  the  ob- 
jective function,  g(X)  is  a  set  of  I  inequality  constraint 
functions,  and  h(X)  is  a  set  of  m  equality  constraint  functions, 

and  X.  and  X.  are  bounds  on  the  design  variables,  referred  to 
11 

as  side  constraints.   The  problem  is  solved  by  creating  a 
single  augmented  Lagrange  function  L(X,X),  where  X  is  a  vector 
of  n  Lagrange  multipliers.   The  optimal  solution  (X*,X*)  is 
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found  by  alternately  solving  a  series  of  unconstrained 
minimizations  followed  by  a  simple  update  of  the  multipliers, 
X. 

In  this  study,  the  unconstrained  minimization  is  performed 
by  the  Davidon-Fletcher-Powell  method  [1]  utilizing  a  combina- 
tion of  the  Golden  Section  and  polynomial  one-dimensional 
search  procedures.   The  Davidon-Fletcher-Powell  method  provides 
good  reliability  and  convergence  properties.   The  Golden 
Section  search  provides  a  reliable  means  of  reducing  the 
search  bounds.   In  most  cases,  polynomial  interpolation 
obtains  a  more  accurate  minimum  than  Golden  Section  for 
sufficiently  narrow  bounds.   The  combination  provided  a 
reliable  and  accurate  one-dimensional  search. 

The  multiplier  method  has  significant  advantages.   As 
will  be  shown  later,  the  inherent  numerical  ill-conditioning 
of  more  common  penalty  function  methods  is  reduced  in  the 
multiplier  method.   Exact  solutions,  not  possible  in  penalty 
methods,  are  also  attainable.   Finally,  multiplier  methods 
are  not  restricted  to  convex  programming  as  are  pure  primal- 
dual  methods. 

There  have  been  numerous  applications  of  the  multiplier 
method  in  the  various  disciplines  of  engineering  design 
optimization.   Imai  [2]  effectively  used  the  method  in 
structural  optimization.   Fax  and  Mills  [3],  while  limited 
to  equality  constraints,  applied  the  method  to  heat  exchanger 
optimization.   Hedderich's  [4]  work  in  heat  exchanger 


optimization  showed  the  need  for  a  method  which  could  treat 
equality  constraints  directly,  in  addition  to  inequality 
constraints.   He  showed  the  Constrained  function  Minimization 
Program,  CONMIN  [5]  to  be  very  inefficient  for  the  equality 
constrained  problem.   It  is  his  work  that  motivated  this 
research  of  the  multiplier  method. 

It  is  the  objective  of  this  research  to  develop  an  opera- 
tional computer  program  applicable  to  the  various  disciplines 
of  engineering  design  optimization.   It  is  not  the  author's 
intent  to  develop  a  universally  superior  optimization  program, 
but  to  develop  a  program  to  be  incorporated  into  a  library 
of  optimization  programs  with  flexibility  in  the  choice  of 
unconstrained  minimization  subprograms,  one  dimensional  search 
subprograms,  derivative  evaluation  subprograms,  and  conver- 
gence criterion. 

The  theory  and  historical  development  of  the  multiplier 
method  is  presented,  concluded  by  a  multiplier  method  algorithm 
for  nonlinear  equality  and  inequality  constrained  problems. 
Experimentation  with  computational  aspects  is  then  summarized. 
Finally,  a  set  of  mathematical  and  engineering  test  cases 
are  solved  demonstrating  the  method's  effectiveness,  followed 
by  conclusions  and  a  discussion  of  results. 


II.   THE  MULTIPLIER  METHOD 

The  development  of  the  multiplier  method  and  a  brief 
mathematical  background  is  presented  in  this  chapter. 
Mathematical  proofs  are  brief  since  computational  aspects 
and  application  of  the  method  are  more  the  concern  of  this 
thesis.   This  chapter  includes  historical  background  of  the 
method,  a  review  of  the  Lagrange  multiplier  and  penalty 
function  methods,  and  a  development  of  the  algorithm  for 
equality  and  inequality  constrained  problems.   Rigorous 
mathematical  development  of  the  method  is  given  by  Hestenes 
[6]. 

A.   BACKGROUND 

The  multiplier  method  was  independently  introduced  in 
1968  by  Hestenes  [7]  and  Powell  [3].  Both  developed  an 
augmented  Lagrange  function  and  solved  a  series  of  uncon- 
strained minimizations  followed  by  a  simple  update  of  the 
multiplier  vector.   Powell  showed  that,  if  the  function  had 
continuous  second  derivatives,  the  method  would  converge 
locally  at  a  linear  rate  while  the  penalty  parameter 
remained  finite.   As  will  be  shown  later,  this  eliminated 
the  inherent  numerical  ill-conditioning  of  classical  penalty 
function  methods. 


Extensive  testing  and  modification  has  been  done  to  the 
method  of  Hestenes  and  Powell.   Miele  et.  al.  [9,10]  applied 
the  method  extensively  to  equality  constrained  problems. 
Tripathi  and  Navandra  [11]  experimented  with  updating  the 
multiplier  after  each  one  dimensional  search  instead  of  after 
each  unconstrained  minimization.   Rockafellar  [12,13]  expanded 
the  method  to  nonlinear  inequality  constrained  problems. 
Global  convergence  of  the  method  was  proved  by  Rockafellar 
[13]  for  convex  programming  problems.   Bertsekas  [14,15] 
provides  convergence  proofs  for  the  general  nonlinear  problem. 

B.   THE  LAGRANGIAN 

The  general  nonlinear  equality  constrained  problem  is 

defined: 

Minimize  f(X)  (2.1a) 

subject  to  h. (X)  =  0  (i=l , . . . ,m<n)  (2.1b) 

The  theorem  associated  with  Lagrange  multiplier  method  states 

that: 

"If  X*  affords  a  local  minimum  to  f (X)  subject  to  the 

constraints  h. (X)=0,  then  there  exists  a  unique  set  of 

multipliers,  A.,  (i-l,....,m)  such  that  if 

m 
L(X,A)  =  f(X)  +  T~       A.h.  (X) ,  (2.2a) 

i-l    1 
then 

m 
VL(X*,A*)  =  Vf(X*)  +  E  A.  Vh. (X*)  =  0  (2.2b) 

i=l   1    1 
and  9 

dZL(X*    A*)    0 

L"(X*fX*)  =  ?1   -  (2-2c) 

dZX± 
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where  V  denotes  the  gradient  of  the  function  and  L" 

denotes  the  second  derivative  of  L"[6]. 
Equations  2.2b  and  2.2c  are  the  necessary   conditions  for 
locally  constrained  minima.   Equation  2.2b  and  the  feasibili- 
ty condition  (Eq.  2.1b)  constitute  the  Kuhn-Tucker  necessary 
conditions  for  optimality  [6].   It  is  assumed  that  f(X)  and 
h. (X)  are  second  order  dif ferentiable  and  that  the  gradients 

Vh. (X)  are  not  zero  at  X*. 
l 

The  problem  can  now  be  stated  in  terms  of  the  equivalent 
classical  Lagrangian. 

Minimize  L(X,A)  (2.3a) 

subject  to  h. (X)  =  0  (i=l, . . . ,m<n)  (2.3b) 

Assuming  the  existence  of  the  saddle  points  of  the  Lagrangian 
L(X,A),  the  condition  exists: 

L(X*,A)  <  L(X*,A*)  <  L(X,  A*)  (2.4) 

The  optimal  pair  (X*,A*)  can  be  obtained  by  first  minimizing 
L(X,A)  respect  to  X,  then  maximizing  L(X  ,  A)  with  respect  to 
A  by  update  equation, 

Ak+1  =  Ak  +  c[h. (Xk) ]  (2.5) 

where  c  is  a  scalar  parameter  (stepsize) ,  k  is  the  iteration 

— k  -k  — k 

number,  and  X   is  the  local  minimum  of  L(X  ,A  ).   The  procedure 

is  repeated  until  convergence  is  attained.   This  is  the  so 

called  "primal-dual"  method. 

Serious  disadvantages  are  encountered  in  the  primal-dual 

method.   First,  the  problem  (Eq.  2.3)  must  have  a  locally 

convex  structure  for  the  dual  problem  to  be  well  defined  and 
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for  Equation  2.5  to  be  meaningful  [16].   Second,  a  large 
number  of  iterations  are  usually  required  to  minimize  L(X,A) 
(Eq.  2.2a)  since  the  ascent  iteration  (Eq.  2.5)  converges  only 
moderately  fast.   Thus,  primal-dual  methods  have  found 
application  in  only  a  limited  class  of  problems  where  minimi- 
zation of  the  Lagrangian  (Eq.  2.2a)  can  be  efficiently  carried 
out  due  to  special  structure,  as  shown  by  Luenberger  [16],  or 
where  the  design  problem  exhibits  a  unique  form,  as  shown  by 
Schmit  and  Fleury  [17]. 

C.   PENALTY  FUNCTION  METHODS 

Penalty  function  methods  have   been  used  extensively  since 

the  mid-1940' s  [18].   They  are  considered  to  be  efficient  for 

inequality  constrained  problems. 

Given   the  nonlinear  inequality  constrained  problem: 

Minimize  f(X)  (2.6a) 

subject  to  g. (X)  £  0  (i=l,...,£)  (2.6b) 

The  general  exterior  penalty  function  is  defined: 

I 
F  (X,c)  =  f  (X)  +  cC  <Mt)  (2.7) 

e  i=l 

where  <j>(t)  is  some  scalar  penalty  function  of  the  constraints 

and  c  is  some  scalar  penalty  parameter.   The  most  common 

2 
penalty  function  is  the  quadratic  <J>(t)=t  /2.   However,  it  may 

be  desirable  at  times  to  use  other  penalty  functions.   In  this 

study  the  quadratic  penalty  function  is  used  such  that 

Equation  2.7  becomes 
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(2.8) 


The  problem  is  now  an  unconstrained  minimization  of  F  (X) . 
The  convergence  of  the  method  is  easily  seen  by  a  simple 
numerical  example. 

Example  1  [2].   Minimize  x  such  that  1-x  <  0. 

Solution:   First,  examine  the  solution  to  the  problem 
by  primal-dual  method.   Equation  2.2a  becomes 

L(x,X)  =  x  +  X(l-x)  (2.9) 

where  A  is  a  real  non-negative  number.   From  the  stationary 
conditions  of  Equation  2.9: 

(x*,A*)  =  (1,1) 
Solving  by  exterior  penalty  method,  Equation  2.8  becomes 

Fp(x,c)  =  x  +  |  p2(x)  (2.10) 

{1-x,    x  >  1 
0   ,    otherwise 
From  the  stationary  conditions,  a  minimum  exists  at 

X'  =  1  -  i 
c 

To  obtain  the  optimal  x*,  a  series  of  unconstrained  minimiza- 
tions are  solved  while  increasing  the  penalty  parameter  c 
toward  infinity.   It  is  apparent  that  as  c->°°,  x'-*x*.   No  exact 
solution  is  obtained  and  the  optimum  is  approached  from  the 
infeasible  region. 
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The  plots  of  F  (x)  vs.  x  and   f(x)  vs.  c  are  given  in  Figs.  1 
and  2,  respectively.   Note  in  Fig.  1,  as  c  gets  larger,  the 
function  changes  more  rapidly  and  the  optimum  becomes  more 
difficult  to  find  regardless  of  the  minimization  technique. 
This  is  a  cause  of  numerical  ill-conditioning  inherent  with 
penalty  methods.   Figure  2  shows  the  asymtotic  convergence 
of  the  method. 

Interior  penalty  functions  have  the  advantage  of 
approaching  the  optimum  from  the  feasible  region  thus 
yielding  a  feasible  solution.   However,  the  penalty  function 
is  discontinuous  at  the  constraint  boundaries.   Also,  the 
same  problems  of  ill-conditioning  and  slow  convergence  exist 
as  seen  by  Example  2 . 

Example  2.   Solve  Example  1  by  the  interior  penalty 

function  method. 

Solution:   One  form  of  the  interior  penalty  function 

is 

F.  ^(x)  =  f(x)  -  -  Y~     7  ;  (2.11) 

mt  c  H — »  g.  (x) 

1=1  3i 

Substituting, 

F.   (X)  =   -  -  -i-  (2.12) 

int         c  1-x 

It  can  be  shown  analytically  that  a  minimum  exists  at 

xl  =  1  +  - 
c 

It  is  again  apparent  that  x'^x*  only  as  c-*°°. 
Figures  3  and  4  illustrate  that  the  same  problems  of  ill- 
conditioning  and  slow  convergence  exist  as  with  the  exterior 
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penalty  method.   An  exact  solution  is  again  not  possible  but 
the  solution  achieved  here  is,  however,  feasible. 

Extended  interior  penalty  methods  [19]  avoid  function 
discontinuity  at  g. (X)  =0  inherent  in  the  interior  penalty 
method.   It  is,  therefore,  less  susceptible  to  ill-conditioning. 
It  does  exhibit  the  same  slow  convergence  of  other  penalty 
methods  due  to  the  requirement  to  increase  c  to  infinity. 
It  is  for  these  reasons  that  the  multiplier  method  is  an 
attractive  alternative. 

D.   THE  EQUALITY  CONSTRAINED  PROBLEM 

The  multiplier  method  can  be  perceived  to  be  a  combined 
primal-dual  and  penalty  function  methods.  Though  they  are 
theoretically  similar,  their  behavior  is  quite  different. 

It  has  been  shown  that  the  original  equality  constrained 
problem  (Eq.  2.1)  is  equivalent  to  the  classical  Lagrangian 
(Eq.  2.3).   Since  Equation  2.3  is  still  an  equality  constrained 
problem,  it  can  be  solved  by  the  usual  exterior  penalty 
function  method.   The  quadratic  penalty  function  is  used  so 
first  derivatives  are  continuous.   Substituting  Equation  2.3 
into  Equation  2.8: 

A(X,A,c)  =  L(X,A)  +  |  \    P?  (x) 

i=l 

m  m,    9 

=  f(X)  +  V       A.h.(X)  +  %     Y2     pf(X)        (2.13) 
4 — i   11       2  u. — '  L  l 
1=1  1=1 
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where  p. (X)  =  h. (X) 
Equation  2.13  is  defined  to  be  the  "Augmented  Lagrange 
function"  for  the  equality  constrained  problem. 

By  the  nature  of  penalty  function  methods,  as  c  goes  to 
infinity,  Equation  2.13  converges  to  the  solution  of  Equation 
2.8.   Concurrently,  by  choosing  an  appropriate  value  of  A 
with  a  suitable  update  formula,  Equation  2.13  can  be  solved 
by  a  series  of  unconstrained  minimizations  to  obtain  a  solu- 
tion to  the  original  problem  (Eq.  2.1). 

The  selection  of  A  can  be  significant  to  the  behavior  of 

the  function  as  seen  in  the  following  two  extreme  cases. 

First,  take  A.=0  for  all  unconstrained  minimizations. 

Equation  2.13  becomes 

m ^ 

A(X,A,c)  =  f(X)  +  %     Li    hi  <x>  (2.14) 

z      i=l 

which  is  the  usual  quadratic  exterior  penalty  function 
(Eq.  2.8).   It  has  previously  been  shown  that  the  function 
only  converges  to  a  minimum  f (X*)  as  c  goes  to  infinity.   It 
has  also  been  shown  to  be  slow  in  converging,  susceptible 
to  numerical  ill-conditioning,  and  to  attain  only  a  near 
optimum  solution  from  the  infeasible  region. 

Next,  consider  the  case  where  A . = A ^ .    At  the  minimum, 
the  stationary  condition  requires  that 
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_  m 

VA(X*,A*,c)  =  Vf(X*)  +  YZ     X-h- (X*) 

i=l 

m 
=  eE  h.(X*)Vh(X*)  =  0  (2.15) 

i=l   x 

The  feasibility  condition  h. (X*)=0  implies  that  Equation  2.13 

is  independent  of  the  value  of  c.   This  leads  to  two  important 

results.   First,  if  the  optimum  A*  is  known  initially,  the 

solution  can  be  obtained  in  one  unconstrained  minimization. 

Second,  since  A(X,A,c)  is  independent  of  c  at  the  optimum,  it 

is  not  necessary  to  sequentially  increase  c  to  infinity  to 

attain  a  solution.   The  second  result  implies  that  a  finite 

c  can  be  chosen,  thus  avoiding  the  inherent  ill-conditioning 

of  the  penalty  methods. 

The  task  is  now  to  find  initial  values  for  c  and  X   with 

sequential  update  formulas  for  each  to  achieve  a  suitable 

rate  of  convergence.   Collecting  terms,  Equation  2.15  becomes 

m 
Vf(X*)  +   F_   [A*+ch.(X*)]  Vh. (X*)  =  0  (2.16) 

1=1 

By  the  Kuhn-Tucker  conditions,  at  the  optimum  (X*,A*). 

m 
Vf  (X*)  +  2_J  x*    7h-  (X*)  =  0  (2.17) 

i=l   1    1 

Thus,  Equation  2.16  reduces  to  Equation  2.17  in  the  limit. 

This  implies  as  update  formula  for  A  such  that 

Ak+1  =  Ak  +  ch. (Xk)  (2.18) 

l       11 

—  k  th 

where  X   is  the  solution  to  the  k    unconstrained  minimization. 

Proposed  initially  by  Hestenes  [6],  Equation  2.13  remains  the 

most  popular  update  formula  for  X. 
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Little  experimentation  has  been  done  with  choosing  and 

— o  — o  —  * 

initial  A  .   As  shown  earlier,  if  A  =A   the  solution  is 

obtained  in  one  unconstrained  minimization.   Obviously,  the 

closer  A   is  to  A*,  the  more  rapid  the  convergence.   A  widely 

accepted  practice  is  to  choose  A?=0  due  to  the  computational 

convenience,  and  because  no  other  multiplier  has  consistently 

proven  more  efficient.   Active  and  violated  constraints  are 

immediately  identified  in  this  case  since  a  non-zero  A.  can 

2  1 

only  occur  for  an  active  or  violated  constraint.   A  constraint 
is  active  if  the  multiplier  is  at  its  optimum  A*^0  when  X-»X*. 
This  eliminates  any  extra  computations  to  check  constraint 
behavior. 

Before  determining  the  choice  of  the  initial  penalty 
parameter  c   and  the  update  formula  for  c,  it  is  first  con- 
venient to  examine  the  convergence  of  the  method  which  is 
directly  related  to  the  choice  of  c.   Various  proofs  of 
linear  convergence  to  a  local  minimum  have  been  developed 
[8,20,21].   Rockafellar  [13]  proved  global  convergence  of  the 
method  for  convex  programming.   Bertsekas  [14,15]  provides 
a  rigorous  proof  of  the  method's  convergence  for  the  general 
nonlinear  problem.   He  compares  the  convergence  of  the 
various  multiplier  methods,  penalty  methods,  and  primal-dual 
methods  [22].   His  results  for  the  Hestenes1  [6]  multiplier 
method  are  summarized  here. 

Recall  the  assumptions  made  for  a  local  minimum  to  exist: 

-k 
Assumption  1.   There  exists  a  local  minimizing  point  X 

of  problem  2.1  which  satisfies  the  second 
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order  sufficiency  conditions  for  an  isolated 

local  minimum,  i.e.  ,  f  and  h.  are  twice 

1 

continuously  dif ferentiable  in  a  neighbor- 

-k  -k 

hood  of  X  ,  the  gradients  Ah. (X  ),  (i=l,...,m) 

are  linearly  independent  and  there  exists  a 

Lagrange  multiplier  vector  A,  such  that 

VL'.X^A^O  and  |  V2L  (Xk,  X X)   |  >0  [22]. 

Assumption  2.   The  penalty  function  cj)(t)  is  twice  continuously 

dif ferentable  in  an  open  interval  containing 

zero  and  cf>"  (0)  >0,  where  c(>"  denotes  the  second 

derivative  [22]. 

2        2    - 

Also,  assuming  the  Hessian  matrices  V  f(X),  V  h. (X),  and  the 

second  derivative  cj)"  are  Lipschitz  continuous,   there  exists 
a  scalar  c*>0,  and  M>0  such  that  for  every  c  >c*  the  function 

If 

L  (X,A)  has  a  unique  minimum  &    (A,c).   Furthermore, 

X"  -X*     -  rAU^ L-LL  (2.19) 

ii  ii         c 

and      [|Xk+1  -  I*M  <  M11^-X*M  (2.20) 

ii  ii        c 

where    A^+1  =  A^  +  c  h± (X^) 

The  notation   | * | I  denotes  the  usual  Euclidean  norm. 

k 

From  Equation  2.20,  if  c  -*-c<°°  for  the  non-decreasing  posi- 
tive penalty  sequence  tc  J  so  as  to  ensure  M/c  ->!,  then  the 

r  — k  \  *  k  j-  • 

sequence  1A  )  converges  to  A   linearly.   If  c   goes  to  infinity, 


A  function  is  said  to  satisfy  a  Lipschitz  condition  of 


order  m  on  the  closed  interval  (a,b)  if  there  exists  a  constant 

-> 

i 
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c  such  that  |f(x  )-f(x  )|<c|x  -x  |  for  all  values  x  ,x   on  (a,b) . 


k  -k 

M/c  ->0,  and  the  sequence  {A  }  converges  superlinearly .   From 

Equation  2.19,  for  sufficiently  large  but  finite  c'  ,  the 

-Jc  -k  k 

sequence  (x  }  converges  to  X*  since  (x  }  converges.   If  c" 

goes  to  infinity,  {X  }  converges  by  the  same  argument. 

-k 
The  rate  of  convergence  of  (X  }  is  represented  by  the  right 

-k 
hand  side  of  Equation  2.19.   Since  i\    }  converges  at  least 

_}- 
linearly,  {X  }  converges  at  least  linearly.   This  has  been 

shown  to  be  significantly  faster  than  pure  penalty  or  primal- 
dual  methods.   Numerical  ill-conditioning  is  avoided  since 
this  linear  convergence  rate  is  achieved  with  a  sufficiently 
large  finite  c. 

The  convergence  criterion  of  Equations  2.19  and  2.20  are 
global  in  nature  since  no  bounds  restriction  has  been  imposed. 
No  restrictions  on  the  convex  or  non-convex  nature  of  the 
problem  are  specified  as  well.   The  global  convergence  of  the 
multiplier  method  is  contingent  upon  the  ability  of  the  un- 

r-k  I 

constrained  minimization  method  to  generate  a  sequence  >-X'  J 

-k  -k 
which  are  well  defined  local  minimums  to  the  function  A(X  ,X  ,c). 

-k  -k 
Naturally,  the  function  A(X  ,A  ,c)  may  have  other  local  mini- 
mums  to  which  the  unconstrained  minimization  method  may  be 
attracted.   Unless  the  unconstrained  minimization  method 
stays  in  the  neighborhood  of  the  same  local  minimum,  the  con- 
vergence argument  is  invalid  and  there  is  no  assurance  that 
the  multiplier  method  will  perform  any  better  or  worse  than 

the  penalty  methods.   It  should  be  noted  that  the  usual 

-  k         th 
practice  is  to  use  the  last  point  X   of  the  k    minimization 
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as  the  starting  point  of  the  (k+l)th  minimization.   This 

— k 
generally  produces  sequences  ix  ;  which  are  close  to  the 

same  local  minimum  of  A(X,X,c). 

The  convergence  rate  of  the  multiplier  iteration  (Eq.  2.18) 

is  linear  with  the  convergence  ratio  essentially  inversely 

proportional  to  the  penalty  parameter  c  [22].   This  fact  is 

strongly  dependent  on  Assumptions  1  and  2.   If  either  of  the 

assumptions  is  relaxed,  the  convergence  rate  may  become 

sublinear  or  superlinear  as  the  following  examples  show. 

Example  3  [22],   Consider  the  scalar   problem  minimize 

2 
x  /2  such  that  x=o  with  an  optimal  solution  x*=0,  A*=0.   In 

this  example  Assumption  2  is  not  satisfied.   For  A^O, 

i  i  3 
4>(t)=|t|  /3/  and  c=l,  Equation  2.13  becomes 


1,.  .,    -1-/(1-4A) 
x  (A,l)  =  ^ 


Starting  at  A  =0,  Equation  2.18  becomes 


Ak+1  =  1-  /(l-4Ak) 


S  k+»( )  =  •*■'  i-e-'  a  sublinear  convergence  rate. 

A^ 

Example  4  [22].   Consider  Example  3  where  cj> (t)=2/| t |/3,  c=l. 

Again,  Assumption  2  is  not  satisfied.   The  solution  to  Equation 

2.13  becomes 
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Starting  at  A  =0,  Equation  2.18  becomes 


xk+l    Ak  +  (-1  +  /l-4Ak) 
It  can  be  shown  that 

lim  J L  =  1 

hence,  a  superlinear  convergence  (order  2). 

Example  5  [22].   Consider  the  problem,  minimize  |x|  /2 

such  that  x=0.   Again  x*=A*=0,  but  Assumption  1  is  not  satis- 

2 
fied.   For  4>(t)=t  /2,  c=l,  the  solution  to  Equation  2.13 

becomes 


xkU/1)  =  tziwpsa 


and  Equation  2.18  becomes 


Ak+1    ^k    (-1  +  /l-4Xk) 


Again, 

lim  I  A    |  _  , 

k*»  (xk)2 

hence,  superlinear  convergence. 

It  can  be  shown  that  the  convergence  rate  is  influenced 

substantially  by  the  rates  of  change  of  the  derivatives 

(curvature)  of  the  primal  function  p(u)  =  min  f(X),  where 

u  =  h(X),  and  the  penalized  primal  function  derivatives 

p  (u) +c  [(J)  (u)  ]  near  u=0.   The  convergence  is  faster  if  the  rate 

of  change  of  Vp(u)  is  small  and  the  rate  of  change  of 

c[V(f)(u)]  is  large  near  u=0. 

24 


In  Example  3,  the  rate  of  change  of  V<J>  (u)  is  small  near 
u=0  and  convergence  is  slow,  while  in  Example  4,  it  is 
large  and  convergence  is  rapid.   In  Example  5,  the  rate  of 
change  of  Vp(u)  is  small  near  u=0,  thus  the  fast  convergence. 

The  following  example  shows  that  in  the  absence  of 
Assumption  1,  Equation  2.18  may  not  lead  to  convergence  for 
any  c>0  when  cj>  is  essentially  quadratic. 


Example  6  [22]  .   Consider  the  problem,  minimize  - 


x 


P 


such  that  x=0  where  l<p<2.   For  any  c>0,  there  exists  a  neigh- 
borhood of  x  =0  such  that  the  augmented  Lagrangian  (Eq.  2.14) 

does  not  have  a  local  minimum  for  any  A  when  <J)  is  essentially 

o  '     2 
quadratic.   This  can  be  corrected  by  using  <£(t)  =  |t|    +  t  /2 

i 

where  l<p  <p. 

The  two  extreme  cases  are  now  examined.   First,  suppose 

_k    - 

"A  -    X* ,    Equation  2.19  becomes 

l|x-k-  X*  ||  <o 

-  k   - 
The  norm  is  always  non-negative;  therefore,  X   =  X*  and  the 

solution  is  reached  in  one  unconstrained  minimization. 

Next,  letting  Ak  =  0  for  all  k,  i.e.,  the  penalty  function 

method,  Equation  2.19  becomes 

i  i  -k   -*.  i  i    M  I  I  A*  |  | 
X   -  X*     <  J—u LL 


-  k  -  k 

In  this  case,  X  +X*  only  if  c  -*°°,  requiring  many  unconstrained 

minimizations  and  a  sublinear  rate  of  convergence. 

It  has  been  shown  that  the  rate  of  convergence  is  directly 

dependent  on  the  penalty  parameter,  c.   The  convergence 

estimates  (Eqs.  2.19,  2.20)  are  valid  for  c  greater  than  some 
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threshold  value,  c*  which  depends  on  (A)  and  the  problem 
data.   In  general,  c*  is  unavailable  and  it  is  impossible  to 
know   apriori  the  range  of  values  of  c  for  which  Equations 
2.19  and  2.20  are  valid  and  imply  fast  convergence. 

A  penalty  parameter  update  sequence  is  required  to  increase 
c  monotonically  with  each  unconstrained  minimization.   As 
c-*-00,  Equations  2.19  and  2.20  will  eventually  become  valid. 
It  should  be  noted  that  large  values  of  c  can  induce  ill- 
conditioning,  making  the  unconstrained  minimization  of  A(X,A,c) 

difficult.   On  the  other  hand,  Equations  2.19  and  2.20  indicate 

-k      - 
faster  convergence  of  {A  }  to  A*  for  large  values  of  c. 

An  update  sequence  is  recommended  whereby  c  is  multiplied 

k  +  i     k 
by  some  constant  y>l,  i.e.,  c    =  yc    .      The  penalty  parameter 

c  is  increased  in  this  manner  to  some  significantly  large  c 

3  max 

Bertsekas  [22]  recommends  y   not  much  larger  than  1  to  avoid 
ill-conditioning  effects  in  the  first  few  unconstrained  mini- 
mization iterations.   This  update  scheme  will  be  used  since 
it  has  been  subject  to  the  most  testing  in  recent  years. 
The  choice  of  an  initial  c  requires  more  experimentation  and 
will  be  discussed  later. 

Other  methods  are  available  for  updating  the  penalty  parameter 
c.   Powell  [7]  suggests  multiplying  c  by  some  constant  (3>1 
only  if  the  violated  constraint,  as  measured  by   | h [X( A , c) ] | | , 
is  not  decreased  by  a  certain  factor  over  the  previous  mini- 
mization, i.e.,  ck+1  =  Sck  if  ! |h[X(Ak,ck)] | |>  Y| |h[X(Xk"1/ck"1] || 
and   c1^"1  =  ck  otherwise,  where  6>1,  and,  y<1  are  some  specified 
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scalars.  This  scheme  generates  a  penalty  parameter  sequence 
that  will  be  constant  after  a  certain  index  and  will  achieve 
convergence  by  virtue  of  enforcing  the  asymtotic  feasibility 
of  the  constraints,  i.e., 

lim  | |h[X(Ak,ck)] ||=0 

Another  similar  update  scheme  is  to  use  a  penalty  parameter 

vector  such  that  c-  is  updated  by  the  Powell  method  only  if 

the  constraint  h. [X(A,c)]  is  violated.   This  case  of  a 

separate  penalty  factor  for  each  constraint  corresponds  to 

merely  scaling  of  the  constraints.   A  simple  modification  to 

Equations  2.19  and  2.2  0  is  used  to  prove  convergence  of 

this  case. 

Finally,  define  a  dual  function  A  (A)  for  the  augmented 

Lagrange  function  A(X,A,c)  such  that 

A  (A)  =  min  A(X,A,c)  (2.21) 

c       x 

k  - 
It  can  be  shown  that  if  X  (a,c)  is  the  solution  to  Equation 

2.21,  then  A  (A)  is  a  twice  continuously  dif ferentiable 
c 

convex  function  with  gradient  given  by 

V,A  (A)  =  Vh[Xk(A,c) ]  (2.22) 

A  C 

where  h[xC(A,c)]  is  the  constraint  vector  (Eq.  2.1b)  evaluated 
at  x*.   Substituting  Equation  2.22,  Equation  2.18  becomes 

Ak+1  =  Ak  +  ck  V,A  (A)  (2.23) 

A  C 

This  relation  shows  that  the  multiplier  iteration  (Eq.  2.18) 
is  an  iteration  of  steepest  ascent  for  finding  the  maximum 
of  the  dual  function,  A  ( A) .   Equation  2.2  3  is  equivalent  to 
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a  steepest  ascent  iteration  for  quadratic  penalty  functions, 

2 
(p(t)    =  t  /2.   The  multiplier  method  can  thus  be  viewed  as  a 

primal-dual  method  with  a  limited  search  for  the  optimal 
Lagrange  multipliers  in  the  dual  space. 

The  Hestenes1  multiplier  method  for  nonlinear  equality 
constrained  problems  will  not  be  summarized.   Quadratic 
penalty  functions  are  used  in  this  algorithm. 

Step  1:   Select  A  =0  and  an  appropriate  penalty  param- 
eter c°>0.   Set  k  =  1. 

-   -k 
Step  2:   Solve  min  A(X,  A  ,c) ,  defined  by  Equation 

-k 
2.13.   Denote  the  solution  X  . 

-k 
Step  3:   Update  X   by  Equation  2.18. 

—k+1    -k  — k  — k 

Step  4:   If  A     =  A  ,  stop.   (X  ,h}    is  the  optimal 

solution.   Otherwise,  go  to  Step  5. 

Step  5:   Set  Ak  =  Ak+1.   If  ck<c    ,  where  c     is 
v  max         max 

some  significantly  large  number,  update  c 

k+1      k 
by  c     =  yc  ,  where  y  is  some  increase 

factor  greater  than  one.   Otherwise,  c  =c    . 

max 

Set  k=k+l.   Go  to  Step  2. 
Note  that  no  assumptions  or  restrictions  have  been  made 
concerning  the  nature  of  the  objective  or  constraint  function. 
If  Assumptions  1  and  2  are  satisfied,  the  solution  will 
converge  to  a  global  minimum;  otherwise,  only  a  local  minimum 
can  be  guaranteed. 
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E.   THE  INEQUALITY  CONSTRAINED  PROBLEM 

The  multiplier  method  has  been  shown  to  be  theoretically 
attractive  for  equality  constrained  minimization  problems. 
In  reality,  many  engineering  problems  involve  inequality 
constraints.   It  is,  therefore,  necessary  to  extend  the 
discussion  of  the  multiplier  method  to  include  inequality 
constraints . 

Consider  the  nonlinear  inequality  constrained  problem: 

Minimize  f(X)  (2.24a) 

such  that  g.(X)  ^  0  (i=l,...,&)  (2.24b) 

where  X  is  the  vector  of  n  design  variables.   Introducing 
slack  variables,  Equation  2.24b  becomes 

g._(X)  +  z?  =  0  (2.25) 

where  z.  is  the  slack  variable  for  the  i    constraint.   The 
problem  is  now  an  equality  constrained  problem  of  the  form 
of  Equation  2.1;  however,  the  number  of  design  variables  has 
increased  to  n+l.      The  augmented  Lagrangian  (Eq.  2.13)  becomes 


a(x,z,a,c)  =  f (x)  +  i2   ^(g^x)   +  XJ) 


i=l 

£— i  ?    2 

+  ?    L^  (gj  (x)  +  XT)  (2.26) 

i=i    x 

If  the  number  of  constraints,  I   is  much  greater  than  the 
number  of  design  variables,  n,  as  is  often  the  case  in 
engineering   design  problems,  the  unconstrained  minimization 

problem  becomes  sizable.   The  scope  of  the  problem  can, 

•  ui      2 
however,  be  reduced  by  eliminating  the  slack  variables,  z^. 
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The  unconstrained  function,  A(X,Z,A,c)  is  first  minimized 
with  respect  to  Z.   For  a  local  minimum  to  exist,  the  station- 
ary conditions 
5A 


9z. 

1 


=  0   (i=l,. . . ,n)  (2.27) 


must  hold.   Differentiating  A(X,Z,A,c) 

H-  =  2Aixi  +  c[gi(X)  +  z2±]     (2z..)  =  0  (2.28) 

A. 
ZiC~c1  +  gi(^)  +  Zi]  =  °  (2.29) 

The  solution  to  Equation  2.2  9  is 

<    -• 

or 

9     A. 

z^  =  — -  -    g. (X) 

1  c    ri 

2 
Since  z.  <  0  is  meaningless,  the  solution  becomes 

2  X- 

zf  =  max  [0,  -gi(X)  -  -±]  (2.30) 

Equation  2.30  shows  that  z.  is  no  longer  an  independent  varia- 
ble.  From  this  equation  it  is  observed  that  if  g. (X)  is  a 

critical  constraint,  z.  =  0.   If  g. (X)  is  non-critical,  z.  > 0, 

l  ^i  l 

Therefore , 

?  A. 

g.  (X)  +  z   =  max  [g±  (X)  ,  ~]  (2.31) 

With  the  slack  variables  eliminated,  the  augmented  Lagrangian 

becomes 

£ 

L 

i=l 


A(X,A,c)    =    f  (X)     +   12      [Xiipi    +    j  ipi] 


X. 

where    4'.    =   max    [g.(X), ]  (2.32) 

1  l  c 
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This  is  referred  to  as  Rockafellar ' s  augmented  Lagrange 
function.   Note  that  A(X,X,c)  has  continuous  first  derivatives 
with  respect  to  c,  but  discontinuous  second  derivatives  at 
gi(X)  =  -X/c.   This  preludes  the  use  of  second  order 
methods,  i.e.,  Newton's  method,  for  unconstrained  minimization 

From  Equation  2.31,  Hestenes'  update  formula  for  X  for 
the  inequality  constrained  problem  becomes 

\l        =   max  [0,  A*  +  c  g± (XK) ]  (2.33) 

The  algorithm  for  Hestenes1  multiplier  method  is  easily  modi- 
fied to  include  inequality  constraints  by  adding  Equations 
2.32  and  2.33  to  steps  2  and  3,  respectively.   Note  that 
since  the  inequality  constrained  problem  is  transformed  to 
an  equivalent  equality  constrained  problem,  the  convergence 
properties  are  identical. 

Example  1  can  now  be  solved  by  the  multiplier  method. 
Example  7  [2].   Minimize  x  such  that  l-x<0. 


Solution:   Equation  2.32  becomes 

f 

A(x, X ,c)  =  x  +  c 


|(l-x)2  +  ^(1-x) ,   if  l-x^0 


x2 

■j— 2  ,      otherwise 

c 
Applying  Hestenes'  algorithm  with  X  =0  and  c  =1,  the  augmented 

Lagrangian  becomes 

A°(x,0,l)  =  x  +  |  J  (1-x)2,    if  l-x>0 

0  ,        otherwise 
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This  has  a  solution  x°=0.   Updating  A  by  Equation  2.32, 
A1  =  max  [0/A°+cg(x'1)  ] 

=  max  [0,  0  +  1(1)]  =  1 
With  A  =1  and  c  =1,  the  next  unconstrained  minimization  of 

A(x,A,c)  becomes 

1  (  2 

A  (x,l,l)  =  x  +    (l-x)V2  +  (1-x)  ,   if  1-x  >  0 

-  1/2  ,   otherwise 

The  solution  becomes  x1  =  1.   Updating  A, 

,2        r«  ,1     1  1  ,  1, n 

a   =  max  [ 0 , A   +  c   g  ( x  ) J 

=  max  [0,  1  +  1(0)]  =  1 

1    0 

Since  A  =A  ,  stop  the  calculation.   The  optimal  solution  is 

(x*,A*)  =  (1,1) . 

Note  the  convergence  of  the  method  even  for  a  constant  penalty 
parameter  c. 

The  functions  A  and  A  are  plotted  in  Fig.  5,  and  the 
objective  value  vs.  c  is  plotted  in  Fig.  6.  Several  facts 
are  evident  from  this  example. 

1.  A  solution  is  obtained  in  a  few  unconstrained  minimi- 
zations. 

2.  Each  unconstrained  problem  is  a  smooth  curve.   Thus 
ill-conditioning  is  avoided. 

3.  Convergence  is  not  asymtotic  and  an  exact  solution  is 
attainable. 

4.  With  an  initial  Lagrange  multiplier  of  zero,  the 
solution  is  obtained  from  the  infeasible  region. 
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Fig.  5.   Augmented  Lagrangian  Function 


exact  objective  value 


OJ 


CO 


Fig.  6.   Convergence  of  Augmented  Lagrangian  Function 
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5.  The  optimal  Lagrange  multiplier  is  obtained. 

6.  Any  starting  point  is  possible. 

The  method,  therefore,  shows  attractive  theoretical  features 
over  pure  penalty  or  primal-dual  methods. 

The  multiplier  method  algorithm  is  now  summarized  for  the 
nonlinear  equality  and  inequality  constrained  problem: 
Minimize  f(X) 
such  that  g. (X)  <  0   (i=l,...,£) 


Algorithm: 


and  h.(X)  =  0   ( j=£+l , ,l+m) 


Step  1.   Choose  an  initial  Lagrange  multiplier  A   =0  and 

an  initial  penalty  parameter  c   >  0. 
Step  2.   Solve  the  unconstrained  minimization  problem 
Min  A(X,A,c)  =  f(X) 

m  « 

L. — '   L  i  r  i    2   l J 
1=1 

m+£  ~ 

+  C      [A.h. (x)  +  f  h;cx)] 

A. 
where  ty.    =  max  [g-(X),  — — ] 

Step  3.   Update  the  Lagrange  multipliers,  A  by 

Xk+1  =  Ak  +  c   max  [g  .  (Xk)  ,  ~]       (i=l ,  .  .  .  ,  I) 

Ak+1  =  A.  +  c  h.(Xk)      (j=m+l,. .. ,m+£) 
]      3      : 

Step  4.   If  Ak  =  Ak+1,  stop.   The  optimal  solution  is 
(X*,A*)  =  (X^,!1^.   Otherwise,  go  to  Step  5. 
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_]^    -k+1 
Step  5.   Set  X      =    X         .   Update  the  penalty  parameter  c  by 

k+1      k 

c     =  yc    . 

If  c     >  c    ,  set  c   =  c    ,  else  set  c   =  c 
max  max 

Set  k  =  k+1;  go  to  Step  2. 

This  is  the  theoretical  algorithm  used  in  the  multiplier 

method.   Minor  variations  are  described  in  the  next  chapter 

which  are  designed  to  improve  the  computational  efficiency, 

reliability,  and  accuracy  of  the  method. 
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III.   COMPUTATIONAL  ASPECTS 

The  effectiveness  of  the  multiplier  method  is  dependent 
upon  the  accuracy  of  the  unconstrained  minimization,  choice 
of  Lagrange  multiplier  update  formula,  penalty  parameter, 
the  form  of  the  penalty  function,  and  the  objective  and 
constraint  functions  themselves.   The  theoretical  effect  of 
the  Lagrange  multiplier  and  penalty  parameters  on  the  method's 
convergence  rate  was  described  in  Chapter  II.   The  require- 
ment that  the  objective  and  the  constraint  functions  satisfy 
Assumptions  1  and  2,  listed  in  Chapter  II. D,  has  also  been 
established. 

This  chapter  will  examine  the  effects  of  the  Lagrange 
multiplier  and  penalty  parameter  from  a  computational  view- 
point.  As  stated  earlier,  the  penalty  function  is  quadratic 
throughout  this  discussion.   The  computational  aspects  of 
the  unconstrained  minimization  program  and  associated  one- 
dimensional  search  subprograms  are  first  examined. 

A.   THE  UNCONSTRAINED  MINIMIZATION  PROBLEM 

A  sequence  of  unconstrained  minimizations  is  required  in 
obtaining  a  solution  by  the  multiplier  method.   Of  the  many 
unconstrained  minimization  techniques  available,  the  Davidon- 
Fletcher-Powell  [DFP]  method  [1]  was  selected  because  of  its 
accuracy  and  rapid  convergence.   Figure  7  compares  the  DFP 
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COMPARISON  OF  UNCONSTRAINED  MINIMIZATION  METHODS 


Banana  Function  Minimization 
F(x)  =  10x1-20x1x2+10x2+x1-2x1+5 


F(x) 


Itr   Nfun 


Steepest  Descent  4.01913  18  262 

Fletcher-Reeves   4.00668  14  211 

CONMIN            4.00294  14  86 

DFP               4.0208  8  99 


(1) 


<s 


in 


ITERATION 
Fig.  7.   Comparison  of  Unconstrained  Minimization  Methods 
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method  to  the  methods  of  steepest  descent,  Fletcher-Reeves 

conjugate  gradient  method  developed  as  part  of  this  study, 

and  the  Fletcher-Reeves  method  contained  in  the  CONMIN 

program  for  the  extremely  ill-conditioned  "banana"  function. 

Note  the  DFP  method  converges  in  fewer  iterations  than  the 

other  methods  with  little  or  no  loss  of  accuracy. 

To  solve  the  general  unconstrained  minimization  problem, 

it  is  necessary  to  first  calculate  a  search  direction,  S, 

initially  the  direction  of  steepest  descent.   A  one-dimensional 

search  is  then  used  to  find  the  minimum  of  the  objective 

function  in  the  search  direction,  S.   The  one-dimensional 

search  computes  a  scalar,  a£  which  minimizes  the  function 

f(X,  ,  +  a,  S,  ) .   The  new  X-vector  associated  with  the  function 

minimum  becomes 

X.  =  X.  ,  +  a,  S, 
k    k-1    k  k 

A  new  search  direction  is  calculated,  and  the  iteration  is 
repeated  until  a  specified  convergence  criteria  is  satisfied. 

The  algorithm  for  the  DFP  method  is  described  as  follows 
for  the  general  unconstrained  minimization  problem: 

Step  1.   Select  an  initial  point  X   and  an  initial  posi- 
tive definite  symmetric  matrix,  H   =  I  (identity 
matrix).   Set  k-1. 
Step  2.   Calculate  the  gradient,  Vf(X  ),  and  set 

S   =  -  H   Vf(X  ) 
o      o     o 

Step  3.   Compute  Xk  =  \_1    +  a£  \ 
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where  a*  minimizes  f (Xj-  +  aS.)  ,  i.e.,  a*  is 
the  solution  to  the  one-dimensional  search. 
Step  4 .   Computer  H    =  H   +  A,  -  B, 

where  AF  =  Vf(X  )  -  Vf(X.  .) 

K  JC  —  -L 

A*  =  \  -  \-i  ■  ~a*h 

-  _    AXAXT 
k    AXTAF 

5     (HAF) (HAF)T 

k  ~  — ^ 

AF   HAF 

Step  5:   Compute  S,+,  =  ~Hv+1vfk/  set  k=k+l,  go  to  Step  3 
The  basic  algorithm  is  modified  to  provide  (1)  scaling  of  the 
variables,  (2)  appropriate  update  and  resets,  and  (3)  a  suita- 
ble termination  criterion. 

A  feature  is  added  which  normalizes  the  search  vector 
such  that 

sk 

Sk      =   - 
NORM    |S,  I 

1  k ' max 

For  excessively  ill-conditioned  problems,  pre-scaling  the 
variables  may  also  be  required. 

Resets  and  updates  of  the  vector,  S  and  the  matrix,  H, 
are  incorporated  to  prevent  method  breakdown  from  round-off 
error  and  other  instabilities.   The  positive  definiteness 

of  H  is  preserved  in  theory  only  if  a*  provides  a  true 

-T    - 
minimum  point,  i.e.,  VF,  ,  S,  =  0.   The  update  of  H  is, 

— T    — 
therefore,  skipped  if  VF,  ,  S,>  £  ,  H^I,  and  the  update 
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was  not  skipped  on  the  previous  iteration.   A  tolerance 
value,  e  =  0.01  was  co-putationally  found  to  provide  a  good 
update  criteria  with  negligible  effect  on  the  convergence 
side. 

Two  checks  were  made  on  S  to  insure  a  valid  search 
direction.   Note  first  that  if  H=I,  the  search  direction  is 
that  of  steepest  descent.   Now  at  iteration  k,  if  S,  •  Vf,  (X)>0, 
the  search  direction  will  not  reduce  the  objective  and  a  new 
search  direction  is  found.   Second,  if  at  the  end  of  the 
one-dimensional  search  f (X,  .) >f (X,) ,  the  function  is  ob- 
viously not  a  minimum  and  a  new  search  direction  is  found, 
In  each  case,  the  H-matrix  is  reset  to  identify,  I,  and  a 
search  in  the  direction  of  steepest  descent  is  performed. 
Finally,  to  maintain  stability  of  the  H-matrix,  it  is  reset 
to  the  identify  matrix  every  NDV+1  iterations,  where  NDV  is 
the  number  of  design  variables. 

The  method  is  terminated  when  the  relative  or  absolute 
difference  of  the  objective  function  is  less  than  0.001  for 
two  consecutive  iterations  or  when  a  preassigned  maximum 
number  of  iterations  is  exceeded.   The  somewhat  strict 
convergence  criterion  was  chosen  since  the  performance  of 
the  multiplier  method  is  highly  dependent  on  an  accurate 
unconstrained  minimization. 

The  DFP  method  has  been  shown  to  be  a  reliable  first 
order  method  with  quadratic  convergence.   It  has  one 
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additional  feature  in  that  it  provides  a  good  approximation 
to  the  Hessian  matrix, 

H  «  32f(X) 


-2 

3X 

Thus,  the  method  provides  many  of  the  advantages  of  a  second 
order,  Newton  type  method,  without  the  tedious  calculation 
of  second  derivatives. 

Finally,  the  DFP  method  requires  subprograms  to  calculate 
gradients  and  perform  the  one-dimensional  search.   Gradient 
calculation  was  done  by  the  first  forward  finite  difference 
method.   While  a  central  difference  method  would  improve  the 
accuracy  of  the  calculations,  the  first  forward  finite 
difference  method  with  a  step  size  of  0.01  has  provided 
acceptable  results  in  experimentation  thus  far.   Discussion 
of  the  one-dimensional  search  subprogram  follows  in  the  next 
section. 

B.   THE  ONE-DIMENSIONAL  SEARCH 

A  combination  of  the  Golden  Section  method  and  polynomial 
interpolation  was  used  for  the  one-dimensional  search.   The 
one-dimensional  search  calculates  a  scalar,  a*  which  minimizes 
the  function,  f(X,  ,  +  a  S,  )  ,   It  has  been  shown  that  the 
accuracy  of  computing  a,  is  critical  to  the  performance  of 
the  unconstrained  minimization  subprogram.   The  Golden  Section 
method  was  used  because  of  its  reliability  and  accuracy.   The 
method  is,  however,  less  efficient  than  polynomial  interpola- 
tion methods  in  many  cases.   The  Golden  Section  method 
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converges  linearly  while  polynomial  methods  have  no  guaranteed 
rate  of  convergence.   It  may  also  require  many  function  evalua- 
tions per  iteration  while  polynomial  methods  require  few. 

The  Golden  Section  method  is  based  on  a  Fibonacci  search 
[16]  which  yields  a  sequence  of  intervals  of  uncertainty  whose 
widths  tend  to  zero  faster  than  any  other  method.   Given  upper 

and  lower  bounds,  x   and  xn,    where  x  denotes  the  scalar,  a. 

u      & 

in  Equation  3.1,  two  interior  points,  x,  and  x~  are  found  such 
that 

Xu  "  X2  =  Xl    Xl 
and  the  ratio 

x2  r- 

—  =  (1  +  /5)/2  =  1.61803 

Xl 
This  is  the  Golden  Section  number.   The  interior  points,  x, 
and  x~  are  defined  as 

X]_  =   (1  -  T)XL  +  TXu 

x2  =  TX£  +  (1  "  T)xu 

where  t  =  (3-/5) /2  =  0.38197 
The  bounds  are  revised  by  a  simple  update  routine  and  one  new 
interior  point  is  calculated  at  each  step.   The  process  is 
continued  until  a  given  relative  convergence  tolerance,  e, 
is  satisfied.   The  algorithm  for  the  Golden  Section  search 
is  given  in  Appendix  A. 

The  relative  convergence  tolerance,  s,    is  defined 
e  =  Ax/(xu  -  x£) 
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Recognizing  that  the  interval  is  reduced  by  the  fraction,  t 
(38  percent)  each  iteration,  the  maximum  number  of  function 
evaluations,  N  per  iteration  is  calculated 

£    =      (1-T) 

or      N  =  In (e) /In (1-T )  +  3 

=  -2.078  ln(e)  +  3 
Note  that  by  defining  e  ,  a  fixed  number  of  function  evalua- 
tions are  performed  each  iteration.   If  the  absolute  conver- 
gence tolerance,  Ax  were  specified,  the  number  of  function 
evaluations  will  vary  each  iteration.   Experiments  have 
shown  an  absolute  convergence  criteria  to  yield  slightly 
better  accuracy  with  nearly  the  same  number  of  function 
evaluations.   Therefore,  an  absolute  convergence  criteria 
was  specified  to  give  the  same  or  better  efficiency  than 
the  relative  criteria  and  with  improved  accuracy. 

The  Golden  Section  method  requires  a  logical  guess  to 

the  first  x   such  that  x   and  x„  bracket  the  minimum, 
u  u       I 

Selection  of  x  was  performed  such  that  if  more  than  two 

u     * 

guesses  were  required,  the  last  three  points  chosen  would 

yield  the  values  x„,  x,  ,  and  x  .   This  provided  a  smaller 
1  V      1'       u        r 

initial  bound  than  if  the  original  x,  was  used  and  also 

saved  one  function  evaluation  since  x,  was  already  found 

for  future  calculations. 

Finally,  a  cubic  approximation  was  performed  using  the 
last  four  points  from  the  Golden  Section  method.   A  better 
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minimum  was  found  in  almost  every  case  at  the  cost  of  only 
one  additional  function  evaluation. 

C.   THE  MULTIPLIER  METHOD 

Experimentation  with  the  multiplier  method  involved  in- 
vestigating the  effects  of  the  penalty  parameter,  c,  the 
update  factor,  y,  and  the  update  formula  for  the  multipliers, 
X.   It  was  desired  to  compute  these  parameters  internally  if 
possible,  such  that  the  method's  performance  would  be  inde- 
pendent of  the  problem.   This  is  necessary  to  provide  a 
method  which  can  be  used  simply  and  reliably  for  many 
engineering  applications.   Effects  of  scaling,  variable 
bounds,  and  perturbed  constraints  were  not  investigated  but 
will  be  discussed  briefly  for  completeness. 

The  choice  of  the  initial  penalty  parameter,  c  ,  can 
have  significant  effect  on  the  efficiency  of  the  multiplier 
method.   No  universally  acceptable  method  has  been  found  to 
select  an  initial  penalty  parameter,  c  .   Two  methods  have 
been  found  to  provide  good  initial  estimates  and  are  recom- 
mended here:   (1)   choose  the  initial  penalty  parameter,  c  , 
of  the  same  order  ot  magnitude  as  the  initial  objective 
function,  or  (2)  choose  the  initial  penalty  parameter  such 

that  |Vf|=c|Vq    I.   The  second  method  requires  the  gradient 
1   '   '  ^max  ' 

of  each  constraint  which  is  not  directly  available  in  the 
current  program.   Both  methods  have  provided  acceptable 
results  for  the  examples  tested.   They  cannot  guarantee  a 
reasonable  convergence  or  a  reliable  answer  in  all  cases. 
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Bertsekas  [22]  recommends  as  increase  factor,  y  not  much 
greater  than  1.   A  y=2  worked  well  in  all  the  cases  tested. 
Reasonable  success  was  achieved  with  y=5,  but  caution  should 
be  used  since  instability  and  ill-conditioning  occurred  when 
c   was  chosen  too  large.   Bertsekas1  recommendation  of  y>l 
to  a  maximum  of  y=2  appears  to  be  adequate  in  most  cases. 

The  maximum  penalty  parameter,  c     is  not  as  critical 

c  J    e  '      max 

to  problem  stability  as  the  initial  c  .   It  was  noted  that 

ill-conditioning  does  occur  for  c     too  large.   A  good 

3  max 

guideline  is  to  choose  c     of  order  five  to  six  times  that 

3  max 

of  c  .   If  ill-conditioning  occurs,  it  will  be  necessary  to 

lower  c    ,  c   or  both, 
max 

Choosing  an  initial  set  of  multipliers,  X   other  than 
zero  was  not  investigated.   The  convenience  of  identifying 
active  and  violated  constraints  without  additional  computa- 
tions when  A.=0  tends  to  overshadow  any  computational  advan- 
tage which  may  be  attained  by  choosing  a  non-zero  initial 
multiplier. 

Normalization  of  the  constraints  is  necessary  in  the 
function  subprogram  to  avoid  breakdown  of  the  multiplier 
method.   Scaling  may  be  required  for  extremely  ill-conditioning 
problems.   Various  scaling  methods  exist,  all  of  which 
provide  some  means  making  the  variables,  objective,  and 
constraints  (or  their  gradients)  equal  to  or  near  the  same 
order  of  magnitude. 
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At  present,  the  algorithm  treats  all  side  constraints 
or  bounds  as  additional  problem  constraints.   Treating  side 
constraints  directly  can  greatly  improve  the  computational 
efficiency,  as  well  as  simplifying  problem  programming. 
Ragsdell  [23]  has  shown  success  in  treating  variable  bounds 
(side  constraints)  directly  but  does  not  offer  any  theoretical 
justification.   Further  research  on  direct  treatment  of  side 
constraints  is  required. 

The  multiplier  method  often  converges  computationally  to 
a  solution  that  is  slightly  infeasible.   This  may  not  be 
desirable  from  the  practical  point  of  view.   Imai  [2]  con- 
sidered the  following  perturbed  problem: 
Minimize  f(X,X) 

such  that  g. (X,X)  £  -  e  (e>0,  i=l,...,m) 
where  £  is  some  small  positive  number.   The  constraints  are 
pushed  slightly  into  the  feasible  region,  so  the  function 
will  terminate  in  the  feasible  region. 

The  augmented  Lagrangian,  L  becomes 
m  9 

l  =  f  +  c  Y2  tVi  + 1  K  (3,2) 

1=1 

X. 
where  \p .    =   max  [g.  +  e,  — —  ] 

Setting  e=0,  Equation  3.2  reduces  to  the  original  Lagrange 

function  for  the  general  inequality  constrained  problem. 

Solutions  obtained  with  e=0  are  only  slightly  infeasible, 

i.e.,  g.-+10~3;  therefore,  accuracy  of  the  original  Lagrange 

function  was  acceptable  even  though  slightly  infeasible. 
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Finally/  a  termination  criteria  was  chosen  such  that 
the  relative  of  absolute  change  in  the  objective  function 
was  less  than  0.001  for  three  consecutive  iterations,  or 
if  a  preassigned  number  of  iterations  was  exceeded.   This 
was  chosen  since  convergence  of  the  objective  function  is 
of  more  interest  than  convergence  of  the  multipliers,  X. 
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IV.   COMPUTATIONAL  EXPERIENCE 

The  multiplier  method  was  used  to  solve  five  test  cases. 
Three  mathematical  examples  were  used  to  demonstrate  the 
method's  computational  performance  as  compared  to  other  methods. 
Two  engineering  examples  were  then  solved  to  show  the  applica- 
tion of  the  multiplier  method  in  engineering  design  optimiza- 
tion.  Results  of  the  multiplier  method  computations  are  given 
in  Table  I.   Table  II  shows  a  comparison  of  the  multiplier 
method  to  other  optimization  methods  for  each  case. 

A.   CASE  1:   THE  CONSTRAINED  ROSEN-SUZUKI  FUNCTION  [5] 

2       2        2        2 

Minimize  f(x)  =  x,-5x..  +x_-5x2+2x3-21x_+x4+7x4+50 

2      2      2      2 
such  that  g,(x)  =  x, +x, +x2~x2+x3+x3+x4-x .-8  =  0 

g2(x)=  x, -X..+2X2+X-  +  2X.-X.-10  -  0 

2       2      2 
g3(x)=  2x]+2x1+x2-x2+x3-x4-5  =  0 

This  problem  was  solved  for  two  individual  cases.   Case  1A 
solves  the  problem  with  constraints  g1(x)  and  g3(x)  as  equality 
constraints.   Case  IB  treats  these  two  constraints  as  in- 
equalities.  The  solution  to  each  case  is  given  in  Table  I. 
This  problem  demonstrates  the  method's  ability  to  solve 
equality  constrained  problems  directly  with  greater  accuracy 
and  efficiency.   As  seen  in  Table  II,  the  method  performs 
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no  better  than  the  exterior  penalty  method  or  CONMIN  for  the 
inequality  constrained  problem.   However,  when  the  equality 
constraints  are  treated  directly,  the  multiplier  method  shows 
significant  improvement.   It  should  be  noted  that  CONMIN  uses 
a  polynomial  interpolation  one-dimensional  search,  which  is 
significantly  more  efficient  than  the  Golden  Section  search 
in  many  cases.   Using  the  polynomial  search  should  make  the 
multiplier  method  more  efficient,  if  reliability  can  be 
preserved.   The  reliability  of  the  Golden  Section  search  made 
it  the  more  desirable  choice  in  the  development  of  the  method, 

B.   CASE  2.   A  SIMPLE  QUADRATIC  FUNCTION  [23] 

2 

Minimize  f(x)  =  4x,-x_-12 

2   2 

such  that  g,(x)  =  25-x,-x2  =  0 

g  (x)=  -10x1+x^-10x2+X2+34  ^  0 


g3  (x)  =  -x1 


< 


0 


g4(x)=  -x2  ^  0 
This  case  was  chosen  as  a  comparison  to  the  Sequential 
Unconstrained  Minimization  Technique,  SUMT  [23],  using  an 
interior  penalty  function  method.   The  problem  solution  is 
given  in  Table  I.   From  Table  II,  it  can  be  seen  that  the 
multiplier  method  is  comparatively  accurate  to  the  SUMT 
method. 
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C.  CASE  3:   THE  PAVIANI  FUNCTION  [2  3] 

.  •  2    2   2 

Minimize  f (x)  =  1000-x  -2x2~x3-x  x  -x  x 

2   2   2 
such  that  g,  (x)  =  x. +x9+x-.-2  5  =  0 

g  (x)  =  8x  +4x  +7x  -56  =  0 

g3(x)  =  -x1    £  0 

g4 (x) =  -x2  <  0 

g5  (x)  =  -x3  <  0 

In  this  case,  it  is  interesting  to  note  that  even  though 
the  accuracy  and  efficiency  is  comparable,  the  penalty  method 
required  a  larger  penalty  parameter,  c  than  the  multiplier 
method.   This  is  a  common  point  and  shows  the  multiplier 
method's  ability  to  avoid  ill-conditioning  by  efficiently 
and  accurately  obtaining  a  solution  with  a  finite  c. 

D.  CASE  4:  THE  THREE  BAR  TRUSS  PROBLEM 

As  a  simple  structural  design  problem,  the  three  bar  truss 

in  Fig.  8  was  considered.   The  problem  was  to  determine  the 

areas  Al,  A2 ,  A3  to  minimize  the  structure  weight,  W.   The 

design  was  subject  to  constraints  -15000<  a.  .<20000  psi  where 

a.   is  the  stress  in  the  truss  member,  i  under  load  condition, 
ID 

j .   An  additional  geometric  constraint  of  A1=A3  was  imposed 
to  maintain  symmetry. 
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THE  THREE  BAR  TRUSS  PROBLEM 


.  3 


MATERIAL:   p  =  0.1  lb/in 

E  =  106  psi 
GEOMETRY:   H  =  10  inches 


LOADS : 


P1  =  P2  =  20,000  lbs. 


STRESS  LIMITS:   -15,000  <  a,   <  20,000  psi 
INITIAL  AREAS:   Al  =  A2  =  A3  =  1.0  sq.in. 


FIGURE  8 . 
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E.  CASE  5:   CANTILEVER  BEAM  PROBLEM 

A  cantilever  beam  with  base,  B,  and  height,  H,  variable 
over  the  length,  L,  was  optimized  to  find  the  minimum  volume. 
The  beam,  as  seen  in  Fig.  9,  was  divided  into  five  sections. 
The  following  constraints  were  imposed:   (1)   bending  stress, 
a   in  each  section  was  not  to  exceed  ±20000  psi.  (2)   deflection, 
5  under  the  load  was  not  to  exceed  ±2.0  inches,  and  (3)  height 
to  beam  ratio,  (H/B)<30.   Additional  side  constraints  were 
also  imposed  on  the  base  and  height.   Separating  the  beam 
into  five  sections  expanded  the  problem  to  one  of  ten  design 
variables  and  37  inequality  constraints.   The  problem  was 
solved  by  the  multiplier  method  with  c=1000  and  y=l.      The 
solution  is  given  in  Table  I,  with  a  comparison  to  the  exterior 
penalty  method  and  CONMIN  in  Table  II.   It  should  be  noted 
that  the  exterior  penalty  method  did  not  obtain  a  solution 
due  to  numerical  ill-conditioning  and  stability  problems. 

F.  SUMMARY 

For  the  cases  involving  equality  constraints,  the  advan- 
tages of  the  multiplier  method  over  penalty  methods  and 
CONMIN  has  been  shown.   The  multiplier  method  has  proven  to 
be  a  suitable  alternative  to  the  other  methods  tested  for 
both  equality  and  inequality  constrained  problems.   While 
the  method  may  have  worked  well  for  the  cases  tested,  it 
cannot  be  guaranteed  to  be  the  most  efficient,  if  even  the 
most  accurate,  in  all  cases.   The  dependence  of  the  multiplier 
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THE  CANTILEVER  BEAM  PROBLEM 


s 


P=10000  lb 


Section  A-A 


rr 


v 


B 


H 


The  beam  is  broken  into  five  segments  of  equal  length 
and  each  segment  is  governed  by  two  design  variables  for 
a  total  of  ten  design  variables. 


Design  requirements 


,6 


MATERIAL:   E  =  30xlCT  psi 

GEOMETRY:   Total  length,  L  =  200  inches 

i  =  1,5 
i  =  1,5 
i  =  1,5 


1  <  H.  <  30  inches 
—   l  — 

0.5  <  B.  <  5  inches 

—   i  — 


H./B.  <  30 
i'  l  — 

STRESS  LIMITS:   (at  the  left  end  of  each  segment) 


a.  =  ~  <  20000  psi 
i    I   — 


i  =  1,5 


FIGURE  9 
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method  on  internal  programming  parameters,  c  and  X,    may 
even  make  the  method  less  attractive  in  some  cases  to  other 
optimization  methods. 
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TABLE  I 
MULTIPLIER  METHOD  RESULTS 

Case  1A:   The  Equality  Constrained  Rosen-Suzuki  Function 
Initial  Obj  =  31.000 

XT   =  (1.0,  1.0,  1.0,  1.0) 
GT   =  (-4.0,  -6.0,  -1.0) 


Final  Obj    =  (6.0075 


XT   =  (0.01607,  1.0285,  1.9799,  -1.018) 
GT   =  (-0.00515,  -0.8896,  0.004387) 


Theoretical  Optimum 

Obj  =  6.0000 

XT   =  (0.0,  1.0,  2.0,  -1.0) 

GT   =  (0.0,  -1.0,  0.0) 
Case  IB:   The  Inequality  Constrained  Rosen-Suzuki  Function 
Initial  Obj  =  31.000 

XT   =  (1.0,  1.0,  1.0,  1.0) 

GT   =  (-4.0,  -0.0,  -1.0) 
Final  Obj    =  6.1178 

XT   =  (0.03206,  1.0334,  1.9509,  -1.054) 

GT   =  (-0.1166,  -0.8921,  -0.1237) 
Theoretical  Optimum  =  Same  as  above 
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TABLE  I  (contd) 

Case  2 :   A  Simple  Quadratic  Function 

Initial  Obj  =  -9.0 

XT   =  (1.0,  1.0) 

GT   =  (23.0,  16.0,  -1.0,  -1.0) 

Final  Obj    =  -31.989 

XT   =  (1.0019,  4.8986) 

GT   =  (-0.0017,  -0.192,  -1.026,  -4.894) 

Theoretical  Optimum 

Obj  =  -32.000 

XT   =  (1.0,  4.8990) 

GT   =  (0.0,  0.0102,  -1.0,  -4.899) 

Case  3:   The  Paviani  Function 

Initial  Obj  =  976.0 

XT   =  (210,  2.0,  2.0) 

GT   =  (-13.0,  2.0,  -2.0,  -2.0,  -2.0) 

Final  Obj    =  961.79 

XT   =  (3.289,  0.2403,  3.7599) 

GT   =  (0.0035,  -0.0048,  -3.286,  -0.241,  -3.761) 

Theoretical  Optimum 

Obj  =  961.715 

XT   =  (3.512,  0.217,  3.552) 

GT   =  (0.0,  0.0,  -3.512,  -0.217,  -3.552) 
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TABLE  I  (contd) 
Case  4:   The  Three  Bar  Truss 
Initial  Obj  =  3.8284  lbs. 

XT   =  (1.0,  1.0,  1.0)  in. 
a    =  14142.2  psi 


a2   =  8284.2 


°31  =  -5858.0 

°12  =  "5858-° 

a22  *  8284'2 

o32  =  14142.2 


Final  Obj    =  2.6  39 


XT   =  (0.7885,  0.4086,  0.7883) 


CT11  =  20001. 


°21  =  14636. 

°31  =  "5364-8 

a21  =  -5364.8 

a22  =  14636. 


°32  =  20005. 


Theoretical  Optimum 
Obj  =  2.632 


XT   =  (0.781,  0.424,  0.781) 

a, ,  and  a__  are  active  constraints 
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TABLE  I  (contd) 
Case  5:   Cantilever  Beam 

Initial  Obj  =  9000  cu.  ft. 

HT   =  (15.0,  15.0,  15.0,  15.0,  15.0) 
BT   =  (3,0,  3.0,  3.0,  3.0,  3.0) 
aT   =  (17778,  14222,  10667,  7111.2,  3555.6) 
(H/B)T  =  (5.0,  5.0,  5.0,  5.0,  5.0) 
6    =  1.0535 
Final  Obj    =  3206.2  cu.  ft. 

HT  =  (26.05,  24.62,  22.70,  19.83,  15.35) 
BT  =  (0.886,  0.796,  0.733,  0.646,  0.506) 
aT  =  (19961,  19904,  19071,  18914,  20111) 
(H/B)T  =  (29.42,  30.92,  30.97,  30.72,  30.32) 
6    =  0.9528 
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V.   CONCLUSIONS 

The  multiplier  method  has  been  shown  to  be  an  accurate 
and  efficient  method  for  solving  problems  in  engineering 
design  optimization.   In  the  particular  cases  tested,  it 
showed  comparable  or  improved  performance  to  other  optimiza- 
tion methods.   In  using  a  finite  penalty  parameter,  c  ,  the 
numerical  ill-conditioning  of  penalty  methods  is  avoided  in 
most  cases.   The  convergence  of  the  multiplier  method  is  at 
least  linear  while  penalty  methods  converge  asymtotically . 
Exact  solutions  are  also  attainable  by  the  multiplier  method. 
The  ability  to  handle  equality  constraints  directly  makes  it 
an  attractive  alternative  to  CONMIN  for  the  equality  con- 
strained problem. 

The  multiplier  method  has  other  advantages  which  make  it 
attractive.   It  can  be  used  as  an  interior  or  exterior  opti- 
mization method,  i.e. ,  the  optimum  can  be  approached  from 
the  feasible  or  infeasible  region.   Any  reasonable  starting 
point  can  be  used.   The  dynamic  selection  of  active  constraints 
when  X°=  0  is  also  a  feature  of  the  multiplier  method. 

Research  of  multiplier  methods  is  far  from  complete.   Its 
use  in  the  various  disciplines  of  engineering  design  are 
endless.   A  multiplier  algorithm  needs  to  be  developed  which 
is  independent  of  internal  computational  parameters.   The 
algorithm  needs  to  be  streamlined  for  easy  application  by  the 
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practicing  engineer.   The  specific  applications  for  which 

the  multiplier  method  is  most  attractive  need  to  be  identified 

The  Augmented  Lagrange  Multiplier  method  is  an  extremely 
useful  program  in  computer-aided  engineering  design.   Its 
applications  today  are  few,  but  its  possibilities  are  endless 
and  invaluable. 
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APPENDIX  A 
ALGORITHM  FOR  GOLDEN  SECTION  ONE-DIMENSIONAL  SEARCH 

The  algorithm  for  the  one-dimensional  search  subprogram 
is  outlined  in  this  appendix.   The  one-dimensional  search  is 
first  performed  by  the  Golden  Section  method.   A  cubic  approx- 
imation is  then  performed  on  the  last  four  values  found  by 
Golden  Section  method.   This  approach  proved  to  be  very 
effective  in  obtaining  an  accurate  solution.   The  one- 
dimensional  search  is  performed  as  follows: 

Step  1.   Specify  the  initial  interval  x   and  x  ,  where 

x  denotes  the  scalar,  a,  in  Equation  3.1. 

Evaluate  the  functions  y.  =  f (x  ) ,  y   =  f(x  ). 

Step  2.   Specify  the  absolute  convergence  tolerance,  Ax. 
Calculate  the  relative  convergence  tolerance, 
e,  and  the  number  of  function  evaluations,  N, 
from  the  equations, 
e  =  Ax/(xu-xj?) 

N  =  -2.078  ln(e)  +  3 

The  value  N  is  calculated  in  floating  point 
arithmetic  and  rounded  off  to  the  next  higher 
integer.   If  N<4,  set  N=4. 
Step  3.   Calculate  x1  and  x2  by  the  equations 

XX  =   (1  "  T)   X£  +  TXu 

X2  =  TX^  +   (1  -  T)   Xu 

where  x  =  (3-/5) /2  =  0.38197 
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Evaluate  y   =  f(x  )  and  y2  =  f(x  ). 
Step  4:   Set  counter,  K=4 .   This  is  because  four  functions 

have  already  been  evaluated.   If  N=4  ,  go  to 

Step  9. 
Step  5:   If  y„  is  greater  than  y  ,  go  to  Step  7. 
Step  6:   y,  is  greater  than  or  equal  to  y~  .   x,  is  a  new 

lower  bound.   Set: 

Xl   =  xl 

^  =  *1 

xl  =  ^2 

*1  =  ^2 

X2  =  TX£  +   d-T)Xu 

y2  =  f(x2) 
K   =  K  +  1 

Go  to  Step  8. 

Step  7:   Y   is  greater  than  y, .   x   is  a  new  upper  bound. 

Set: 

Xu  =  X2 

^u  =  y2 
X2  =  Xl 

y2  =  y1 

xl  =  (1"T)X£  +  TXu 
y1   =  f(xL) 

K   =  K  +  1 
Step  8:   Check  convergence.   If  K>N,  go  to  Step  9;  else 

go  to  Step  5. 
Step  9:   Do  cubic  approximation  with  values  x.,    x^,    x2   x^. 

^0   =    f(x£'  V  X2'    Xu] 
Step  10:  Pick  best  of  y£,  y1 ,  y2,  yc- 

Y  =  min  (y£,  y±,    Y  2'    Yc) 

X  =  corresponding  x«f  x^,    x2  or  xc 

Y  is  the  optimum  function  value  at  location  X.  Stop 
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